Inter‐ and intraspecific variation in mycotoxin tolerance: A study of four Drosophila species

Abstract Many mycophagous Drosophila species have adapted to tolerate high concentrations of mycotoxins, an ability not reported in any other eukaryotes. Although an association between mycophagy and mycotoxin tolerance has been established in many Drosophila species, the genetic mechanisms of the tolerance are unknown. This study presents the inter‐ and intraspecific variation in the mycotoxin tolerance trait. We studied the mycotoxin tolerance in four Drosophila species from four separate clades within the immigrans‐tripunctata radiation from two distinct locations. The effect of mycotoxin treatment on 20 isofemale lines per species was studied using seven gross phenotypes: survival to pupation, survival to eclosion, development time to pupation and eclosion, thorax length, fecundity, and longevity. We observed interspecific variation among four species, with D. falleni being the most tolerant, followed by D. recens, D. neotestacea, and D. tripunctata, in that order. The results also revealed geographical variation and intraspecific genetic variation in mycotoxin tolerance. This report provides the foundation for further delineating the genetic mechanisms of the mycotoxin tolerance trait.


| INTRODUC TI ON
Many species within the genus Drosophila have radiated to use a wide variety of hosts for feeding and breeding (Markow, 2019).
These hosts are chemically and phenologically distinct and include fruits, flowers, cacti, slime fluxes, and mushrooms. These adaptations involve genetic and genomic changes. One such adaptive radiation is mycophagy. Many Drosophila species are mycophagous, and mushrooms appear to provide all the essential components of an insect diet (Courtney et al., 1990). However, some mushroom species also contain highly lethal compounds (mycotoxins) to protect themselves from mycophagy (Stump et al., 2011). Although the toxic mushroom species are fewer in number as compared to the nontoxic mushroom species (Graeme, 2014) and only constitute a small portion of the potential diet, many mycophagous Drosophila species can tolerate high concentrations of α-amanitin, which is the most potent mycotoxin (Jaenike et al., 1983;Lacy, 1984;Spicer & Jaenike, 1996;Stump et al., 2011).
Seventeen mycophagous Drosophila species from five species groups within the immigrans-tripunctata radiation have been shown to tolerate mycotoxins (Izumitani et al., 2016;Scott Chialvo & Werner, 2018). These species groups are tripunctata, testacea, cardini, bizonata and quinaria. Very little is known about the feeding habits for the species groups cardini and bizonata. The tripunctata species group comprises 83 species (O'Grady & DeSalle, 2018), and for most species in this group, larval feeding substrates have not yet been determined. The testacea species group contains four species, all of which are mycophagous, whereas 34 species belong to the quinaria group, most of which are mycophagous (O'Grady & DeSalle, 2018;Scott Chialvo et al., 2019). The quinaria species group is of particular interest as mycophagy has been gained and lost multiple times within this group. Furthermore, the loss of mycophagy has been followed by a loss of toxin tolerance without an evolutionary lag (Spicer & Jaenike, 1996), suggesting that mycotoxin tolerance is probably a costly trait.
Although the association between mycophagy and mycotoxin tolerance in certain Drosophila species was established almost three decades ago, the genetic mechanisms involved in the tolerance are mainly unknown. The most lethal mycotoxin, found in the notoriously deadly Amanita mushrooms, alpha-amanitin, binds to RNA-polymerase II (RNA-pol II) and hinders its function. Jaenike et al. (1983) observed that the tolerance mechanism apparently did not involve target modification of RNA-pol II. Another study demonstrated that Phase I detoxification enzymes (Cytochrome P450s) might be conferring mycotoxin tolerance in some but not all mycophagous species (Stump et al., 2011). Apart from these few reports, the understanding of the genetic basis of mycotoxin tolerance has remained inadequate.
To identify mechanisms that confer mycotoxin tolerance, we must understand the genetic architecture of the trait. To achieve this goal, we consider the following questions: (1) Does the mycotoxin tolerance trait show intraspecific genetic variation? (2) Do different species demonstrate variation in the extent of mycotoxin tolerance?
To address these questions, we have performed mycotoxin tolerance assays on multiple isofemale lines of four species within the immigrans-tripunctata radiation. Figure 1 provides images of the four species: D. falleni, D. recens, D. neotestacea, and D. tripunctata (Werner et al., 2020). Figure 2 shows each species in the phylogenetic context. Drosophila tripunctata belongs to the tripunctata species group (Clade A), D. neotestecea belongs to the testacea species group (Clade B), D. recens and D. falleni belong to the quinaria species group, Clade C1, and Clade C2, respectively. Each of these four species represents a major clade of the immigrans-tripunctata radiation and is known to be mycotoxin tolerant.
Previous research on mycotoxin tolerance in Drosophila has focused on α-amanitin, the toxin found at high concentrations in some Amanita mushroom species (Garcia et al., 2015;Jaenike et al., 1983;Jaenike, 1985;Lacy, 1984;Spicer & Jaenike, 1996;Tuno et al., 2007;Stump et al., 2011). However, toxic mushrooms contain a myriad of different toxins (Yin et al., 2019). Therefore, studies based on a single toxin in isolation have a drawback. They cannot account for the potential synergistic and antagonistic interactions among different mycotoxins found in wild toxic mushrooms. In this study, we have used a natural-toxin mix (Scott Chialvo et al., 2020) as a source of mycotoxins. First, using this natural-toxin mix extracted from Amanita phalloides mushrooms ensures a mycotoxin representation similar to that found in the wild, as it contains α-amanitin, β-amanitin, ɤ-amanitin, amanin, amanullin, phallacidin, phallisacin, phalloin, phallisin, and phalloidin (Scott Chialvo et al., 2020). The second reason is that α-amanitin is expensive, and therefore, the use of the natural-toxin mix proved to be cost-effective for this large-scale study. for fly collection within each location spanning over 3-5 square kilometers. The species and the sex of the captured flies were identified, and isofemale lines were set up by adding one wild-caught female with one wild-caught male from the same species and location and collecting their progeny (David et al., 2005). The established isofemale lines were maintained on a diet of Carolina Biological Formula 4-24 Instant Drosophila Medium supplemented with finely ground, freeze-dried Agaricus bisporus mushrooms (Oregon mushrooms, OR) at a ratio of 33.28:1 w/w, and a dental roll was added to the food vial as a pupation site. The standard conditions for maintenance and experiments were 22°C and a 14 h:10 h (L:D) photoperiod at 60% humidity. The authors note here that the isofemale lines were maintained in the laboratory for at least over a year before the experiments were conducted. F I G U R E 1 Representative images of the four organisms used in this study

| Mycotoxin tolerance assays
Basic food was prepared by mixing 28.3 g freeze-dried A. bisporus mushrooms (Oregon mushrooms, Oregon) with 941.9 g Carolina 4-24 Instant Drosophila Medium and grinding them together into a fine powder. For mycotoxin tolerance assays, clean glass vials were filled with 250 mg of basic food.
The natural-toxin mix was provided by Dr. Clare Scott-Chialvo (Scott Chialvo et al., 2020), which contained methanol as eluent. To account for this methanol, one milliliter of 0.56% methanol solution was added to the control vials containing 250 mg of basic food. The mycotoxin vials were prepared by adding 1 ml of the natural-toxin mix (100 μg/ml of known amatoxins) to the vials containing 250 mg of basic food. Both control and mycotoxin vials were weighed and subjected to vacuum evaporation for 96 h at room temperature to remove methanol from the vials. The loss in weight (in grams) in the vials was replenished with the appropriate amount (in ml) of sterile distilled water. The optimal duration of vacuum evaporation was identified using preliminary studies and 96 h of vacuum evaporation showed survivorship that was comparable to vials without methanol.
Water agar plates were prepared using 15 g Bacto Agar (Sigma Aldrich) in 500 ml of distilled water and adding Tegosept to a final concentration of 0.1% and poured into 30 mm Petri-plates.
These plates snugly fit the plastic bottles that were used to make egg-laying chambers. Tiny holes were punched into these plastic bottles for aeration. Equal amounts of dry yeast and freeze-dried mushroom powder were mixed together with autoclaved distilled water to prepare a paste (prepared fresh daily). A drop of this paste was applied to the water agar plate. Recently eclosed males and females of each isofemale line were transferred to egglay chambers and allowed to oviposit at 22°C and a 14 h:10 h (L:D) photoperiod at 60% humidity. The next day, the plates were replaced with fresh plates, and the water agar plates with oviposited eggs were allowed to hatch at 22°C and a 14 h:10 h (L:D) photoperiod at 60% humidity. The hatched first-instar larvae were used for the experiments. Pilot studies were performed to identify the optimal larval density for each species. As a result, 15 first-instar larvae were added to each vial in the case of D. falleni, D. recens, and D. tripunctata, whereas 20 first-instar larvae were added to each vial for D. neotestacea. The experiments were conducted in triplicates. Each experiment was conducted on consecutive days to generate three replicates for each of the 10 isofemale lines/ location/species for two treatments (control and mycotoxin).

| Development time, thorax length measurements, and survival
The vials were checked daily to record the time to pupation, survival to pupation, time to eclosion, and survival to eclosion. The eclosed flies were collected within 24 h by light CO 2 anesthesia, sexed, and placed laterally to measure the thorax length. The thorax's anterior margin length to the scutellum's posterior tip was measured and recorded as the thorax length. The thorax length of the eclosed flies was measured to the nearest 0.025 mm with an Olympus SZX16 dissection microscope fitted with an Olympus DP72 camera, using the ImageScan software (Hasson et al., 1992).
The eclosed females were used for the fecundity assays, and the eclosed males were used for the longevity assays. The experiments were terminated after ensuring that no new flies had emerged for four consecutive days.

| Fecundity assays
Only female flies were used for the fecundity assays. Female flies eclosed from the mycotoxin tolerance assay vials were labeled appropriately and maintained individually in food vials for 3 days as virgins. They were then transferred individually into a fresh food vial with three 3-day-old virgin males from the laboratory stocks of the same isofemale line. These parent flies were transferred to a new vial every 3 days. After 15 days, the adult flies were removed. The offspring of the females that survived the full 15 days were counted to provide an estimate of fecundity. We note that this assay cannot be used to evaluate egg-to-adult survival (Dyer & Hall, 2019).

| Longevity assays
Only male flies were used for the longevity study. Male flies eclosed from the mycotoxin tolerance assay vials were maintained individually in tiny 5-ml glass vials containing approximately 250 mg of the basic food used to create the mycotoxin tolerance assay vials. The vials were checked every alternate day to record any dead flies, and the remaining flies were transferred to fresh food vials every 2-3 days.

| Statistical analyses
All statistical analyses were done using R version 3.6.1 (https://www.r-proje ct.org/found ation/) and R Studio version 2021.09.2 + 382 (https://www.rstud io.com/). We used the linear mixed model (LMM) and the generalized linear-mixed effects model (GLMM), implemented in R package 'lme4' (Bates et al., 2015), to determine the independent variables that can explain the variation in survival, development time, body size, fecundity, and longevity.
We modeled pupal and survival to eclosion using a binomial linear mixed model with the logistic link function. We used the linear mixed model to model development time, fecundity, longevity, and thorax lengths. The development time and thorax lengths were analyzed for each sex separately. To detect whether mycotoxin tolerance shows interspecific variation, we first fitted a GLMM that includes the main effects, the two-way interactions, and the three-way interaction of the species, the treatment, and the location as the fixed effects. The likelihood test (LRT) was used to test if the three interaction was significant. The final model includes the main effects and the twoway interactions of the species, the treatment, and the three-way interaction only if the three-way interaction is significant (in other words, the p-value from the LRT for the three-way interaction is less than 0.05). In all models, the isofemale lines and the replicate vials were included as the random effects. To check the sufficiency of the model, the scatter plots of the deviance residuals against the predicted values were generated, and the dispersion parameter was estimated based on the ratio of the sum of squared deviance residuals and the degrees of freedom of the model if the binomial linear mixed model was used.
To evaluate the effect of toxin treatment, we fitted either a binomial linear mixed or linear mixed model to assess whether the treatment affects the seven gross phenotypes. In all models, the main effect of the treatment was the only fixed effect, and the isofemale lines and the replicate vials were included as the random effects. The analysis was conducted for each of four species and seven gross phenotypes separately. Among seven gross phenotypes, the development time of eclosion and the thorax length of eclosion were analyzed for the males and the females separately. Therefore, 36 p-values were obtained. To account for the multiple testing adjustment, both the p-value adjusted using the Bonferroni correction and the false discovery rate (FDR) calculated with the Benjamini-Hochberg method (Benjamini & Hochberg, 1995) were presented.
The FDR < 0.05 was used as a cutoff for significant results.
To identify the extent of tolerance for each isofemale line, the binomial linear mixed model was performed on each isofemale line with the replicate vial as a random effect, and the lines were segregated based on their p-values. High-tolerance lines were identified as those in which no significant difference in survival between the control and the mycotoxin treatments was observed or where the survival was significantly higher in the mycotoxin treatment.
Isofemale lines with significantly low survivorship on the mycotoxin treatment (p-value <.05) were categorized as low-tolerance lines.
For the scope of this study, high tolerance is defined as the ability of an isofemale line to survive in the presence of the natural-toxin mix (100 μg/ml of known amatoxins).
For intraspecific variation, before model fitting, we pruned the data to exclude isofemale lines where only one data point was observed per treatment. This exclusion allowed us to eliminate data that could not estimate variation within an isofemale line. We as-

| Interspecific variation
For survival to pupation and survival to eclosion, the p-values for the three-way interaction were 0.273 and 0.431, respectively. Either of them achieved the significance level of 0.05. Therefore, the models with all the main effects and the two-way interactions of the species, the treatment, and the location were used. Tables 1 and 2 present the p-values from the LRT for each term in the models for survival to pupation and survival to eclosion, respectively. The least square estimates of estimate and bounds (95% confidence intervals) of the logarithmic of the odds ratio between each pair of species stratified by the treatment are presented in Tables S1 and S2, respectively. We observed a significant interspecific variation in mycotoxin tolerance for survival to pupation (p-value <.001, Table 1) and The significant treatment effects are presented in Table 3. The results for all four species and seven gross phenotypes can be found in Table S3. Interestingly, the effect of mycotoxin treatment followed a similar trend (Table 3)

| Intraspecific genetic variation
In all the four species, the traits that were significantly affected by the interaction between treatment and the isofemale line are tabulated in Table 4. All four species (D. falleni, D. recens, D. neotestacea, and D. tripunctata) showed intraspecific genetic variation in mycotoxin tolerance for survival to pupation and pupal development time.
Additionally, each species showed intraspecific genetic variation in other traits as well (

| Geographical variation
The location and treatment interaction significantly affected survival to pupation (p-value <.05) and survival to eclosion (p-value <.05) (Tables 1 and 2). Tables S45 and S46 include the estimates, bounds, and effect sizes. Generally, isofemale lines from the GSM location showed a poorer survival due to mycotoxin treatment as compared to the ESC location.
While performing statistical analyses for intraspecific genetic variation, we also observed that the location and treatment interaction significantly affected the development time in D. tripunctata.
Males and females of D. tripunctata isofemale lines from ESC showed a significant increase in development time as compared to their GSM counterparts (false discovery rate (FDR) < 0.05, each; Figure 11, Table 4, and Table S8).

| DISCUSS ION
This study provides three key findings. Firstly, it demonstrates sig- Within the quinaria species group, D. falleni and D. recens split ~20 million years ago (Mya) (Izumitani et al., 2016). In our study, D. falleni appears to be more mycotoxin-tolerant than D. recens. A comparison between D. falleni and D. recens has been made previously (Jaenike, 1985), where it was observed that D. falleni larvae could not survive to adulthood at α-amanitin concentrations above 500 μg/ml of   (Izumitani et al., 2016). We suggest that different selective pressures that have brought about the divergent evolution and speciation among these four species have also affected mycotoxin tolerance. It would be interesting to identify the genetic and genomic changes that may have altered the extent of mycotoxin tolerance among these species.
We used the isofemale line technique to investigate the intraspecific genetic variation of mushroom toxin tolerance (Hoffmann & Parsons, 1988). This technique is based on a simple concept that when isofemale lines from wild-collected females are established, and their progeny is maintained under similar laboratory conditions, the variation observed among the isofemale lines is primarily genetic.
Phenotypic variation for a trait among these genetically distinct isofemale lines would indicate intraspecific genetic variation attributed mainly to segregating alleles at multiple loci (Mackay, 2010).
Our study found phenotypic variation in mycotoxin tolerance among isofemale lines in all four species, providing evidence for intraspecific genetic variation for mycotoxin tolerance. Intraspecific genetic variation in mycotoxin tolerance has been reported previously in D.
tripunctata (Jaenike, 1989). Our results confirm their findings and expand the dataset to include three additional species (D. falleni, D. recens, and D. neotestacea). This study provides the groundwork for further studies to calculate the heritability and identify the genetic architecture of the mycotoxin tolerance trait.  (Zhang et al., 2018) to the chemical composition of phenolics in plants . In essence, if a trait shows differences between populations from different geographical TA B L E 4 The traits that show the significant effect of interaction between treatment and Isofemale line (false discovery rate < 0.05) and the corresponding chi-square statistic and its degrees of freedom, the p-value without the adjustment, the p-value after the Bonferroni correction, and the false discovery rate (FDR) calculated based on the Benjamini and Hochberg's method that play a critical role in the evolution of the mycotoxin tolerance trait.
All aforementioned conclusions are based on the results from the linear mixed model or the binomial linear mixed effects model with the logistic link function implemented in R package 'lme4' (Bates et al., 2015). In all models, the replicate vials and/ or the isofemale lines were included as the random effects. It is well known that the binomial linear mixed model may not be sufficient due to overdispersion. Among 19 estimated overdispersion parameters, 16 of them were between 0.93 and 1.52, and only three of them were greater than 1.60 (1.61, 1.86, 2.02, see Figure   S34), indicating that the overdispersion was not a serious problem. The scatter plots of the deviance residuals and the predicted values (see Figures S1-S33) show that most of the model fittings were adequate.
In conclusion, our study identifies interspecific and intraspecific variation in mycotoxin tolerance and demonstrates geographical variation in mycotoxin tolerance. Formal analysis (supporting); software (equal); validation (supporting); visualization (supporting); writing -review and editing (supporting). Thomas Werner: Conceptualization (equal); data curation (supporting); formal analysis (supporting); funding acquisition (lead); investigation (supporting); methodology (equal); project administration (lead); resources (equal); supervision (lead); validation (supporting); visualization (supporting); writing -original draft (supporting); writing -review and editing (lead). and Dr. Laura Reed for guiding the data analyses. The authors also thank Patrick McKee for his assistance in the data collection. We also thank the three anonymous reviewers who provided excellent comments and suggestions that have improved the manuscript.

CO N FLI C T O F I NTE R E S T
There are no competing interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
Phenotype data and the code for running analyses are available in Dryad, https://doi.org/10.5061/dryad.2ngf1vhr7.